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SUMMARY 


This report is concerned with the development of a general method for calcu- 
lating unsteady heat transfer on turbine blades. It Is based on the numerical 
solution of the boundary-layer equations for laminar, transitional and turbu- 
lent flows. A novel procedure has been developed to account for the movement 
of the stagnation point caused by blade-passing wakes and has been applied to 
the stagnation region of three model flows with results which confirm Its 
validity for laminar flows. It has also been applied to an experimental 
arrangement In which the average Nusselt number has been reported for turbulent 
flows and the results show considerable promise. 
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1.0 INTRODUCTION 


Gas turbine designers are constantly expanding the pressure and temperature 
operating thresholds of modern turbines In order to achieve greater power 
densities. This In turn requires careful thermal-structural analysis In order 
to ensure that these turbines, which deliver such excellent performance, will 
also run for a long time before maintenance Is required. The costs associated 
with hot section maintenance account for approximately 60 percent of total air- 
craft engine maintenance [1] and the turbine Industry Is aggressively pursuing 
programs to Improve the predictability of turbine life [2]. In addition, NASA 
has a durability technology program, known as the Hot Section Technology (HOST) 
project [3,4] directed at Improving the predictive tools for analyzing aircraft 
gas turbines and another program directed at Improvement of the durability of 
the Space Shuttle Main Engine (SSME) [5]. 

The ability to predict the metal temperature distribution Is a key factor In 
life prediction. Coincident with the beginning of the HOST project, Stepka [6] 
estimated that the ability to predict steady-state metal temperature In an 
operating engine was within 1 00K , and that this could be refined to 50K by 
testing prototypes. This level of thermal accuracy Is estimated to lead to an 
order of magnitude uncertainty In life prediction [7]. 

Periodic reviews of gas turbine heat transfer [8-10] suggest that progress Is 
being made In measuring and analyzing the complex flows associated with steady- 
state heat transfer In gas turbines. Research on unsteady heat transfer In 
turbine passages Is, however, just beginning to appear. It Is premature to 
evaluate the effects of wake-generated unsteady heat transfer on turbine dura- 
bility but the research data to provide the necessary thermomechanical loads 
are emerging rapidly [11-20] with most of the Information of aerodynamic prop- 
erties acquired In large low-speed turbines. Drlng, et al. [11] explored the 
nature of the boundary-layer response to turbine wakes and the unsteady pres- 
sure loading using thln-fllm surface sensors. Hodson [13] and Binder, et al. 
[14] mapped the convection of wakes through the rotor and flow properties In a 
rotor at 8850 rpm have also been reported In [14]. 

The nature of the blade-passing problem can be explained In relation to Fig. 
la which shows that the rotor blade passes through the wake of the upstream 
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(a) 


(b) 


Figure 1. Flow configuration, (a) Turbine stage, (b) Simulation. 

stator blade and generates a wake which affects the onset flow to the second 
stator. In effect, the blades of the rotor and second stator are subjected to 
a freestream velocity which varies periodically with time and with a random 
turbulence fluctuation superimposed. These effects may be simulated by the 
arrangement of Fig. lb which shows one cylinder of a number arranged on a wheel 
which rotates upstream of a turbine blade. 

Most relevant heat-transfer data have been obtained In short-duration facil- 
ities [15-18] or by using short duration techniques [19,20]. Thln-fllm 
resistance thermometers provide Information of surface temperature and are 
combined with a one-dimensional semi-infinite heat-conduction analysis to 
yield time-resolved heat flux. Dunn and coworkers have pioneered this tech- 
nique and, In their most recent work [15,16], have begun to obtain time- 
resolved heat flux data in a real turbine stage at fully-scaled conditions. 

At the University of Oxford [17,18], rotor wakes are simulated with a spoked 
wheel rotating In front of a linear cascade and the results show that the wakes 
Impinging on the airfoil can cause the local boundary layer to undergo transi- 
tion and reverse transition at the frequency of the passing wakes. O'Brien et 
al. [19] used the transient technique In a steady-running rig, also with a 
spoked wheel wake generator, and showed that the time-average of the time- 
resolved heat transfer agreed closely with conventional steady-state data. 
Morehouse and Slmoneau [20], used the same rig, see Fig. 2, and showed that It 
was necessary to separate the heat transfer effects due to the unsteady dis- 
turbances of periodic wakes from those generated by conventional turbulence. 

The recent advances In computational fluid dynamics have made It possible to 
predict blade pressure distributions, detailed boundary-layer and wake 
profiles, skin-friction and heat-transfer coefficients by solving the steady, 
full Navler-Stokes equations, see for example Shamroth [21]. Also Ral 
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Figure 2. Experimental set-up of O'Brien et al. [19]. 

[22] has attempted to reproduce Drlng's measuremments by constructing a Navler- 
Stokes computer code with a complex arrangement of patched and overlaid grids. 
Navler-Stokes solutions have the advantage of enabling the prediction of both 
attached and separated viscous flows and hence can be expected to become a 
valuable engineering method as soon as computation times can be reduced and 
questions concerning the computational accuracy can be resolved. 

In this report an alternative approach Is pursued which has proven reliable 
and powerful In many external flows. It Is based on the use of Interactive 
boundary-layer theory which requires the development of viscous and Invlscld 
methods and their coupling by special techniques such as those described In 
[23-26]. The Invlscld method must be able to cope with the periodic onset 
velocities of Fig. lb and their effects on the blade flow from the stagnation 
region to wake. The boundary-layer method for the resulting unsteady flows has 
to account for the reversals associated with the movement of the stagnation 
point with time and space as well as flow reversals and separation that may 
occur In the downstream region. The problem of flow reversal can be solved 
without Interaction between Invlscld and viscous flow solutions but a special 
numerical procedure Is required to deal with the changing flow direction. 

The work described here Is directed towards the development of a general method 
for calculating unsteady heat transfer on turbine blades. As an essential 
preliminary we describe a general boundary layer method for calculating heat 
transfer for prescribed freestream conditions and test It for three model 
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problems which emphasize the stagnation region. We assume that the external 
velocity distribution Is represented by a function 

u ( x , t) /u = A(*,t)[5 - B( t) ] (1) 

which allows the variation of the stagnation point and the freestream velocity. 

The first model problem corresponds to flow over a thin ellipse In a constant 

freestream velocity field with a uniform wall temperature and with a local 

external velocity which changes with angle of attack. The second model problem 

corresponds to an Initially steady flow which accelerates from u /u = l at 

6 00 

t = 0 to u /u = 2(1 - 1) at large times with corresponding movement of the 
e °° 

stagnation point from 0 to 1 . The third model problem corresponds to the flow 
of [20], see Fig. 2, In which the flow over a circular cylinder with a uniform 
wall temperature Is subjected to an external velocity distribution which 
changes with time due to the Interference of a series of cylinders which are 
attached to a rotating wheel. Measurements of local flow properties and of 
wall heat transfer have been reported In [20] and allow quantitative assessment 
of the calculation method. 

The basic equations, Initial and boundary conditions are considered In the 
following section which Is followed by a description of the numerical method 
used to solve the boundary-layer equations. Section 4 presents the results and 
discussion and the report ends with a summary of the more Important 
conclusions. 


| 
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2.0 BASIC EQUATIONS 


For two-dimensional, incompressible time-dependent laminar and turbulent 
flows, the boundary-layer equations and their boundary conditions are well 
known. For specified wall temperature, no mass transfer and with eddy viscos- 
ity, c , and turbulent Prandtl number, Pr. , concepts they can be written 
m t 

in the form [27]: 



3u 3v . 

— + T~ a 0 

3x 3y 


(2) 

du 3u 

T7 + U 7 — + 

at ax 

, 3u 3u . 

3u e e 3 

v = -rr~ + u ^ — + v 

3y at e 3x 3y 

th 

<b 1 9y> 

(3) 


3T 3T v 3 , h 3T. 

3x v 3y = Pr 3y l °2 ay' 


(4) 

ll 

O 

u = v = 0, T - yx) 


(5a) 

y = 4, 

u = u ( x , t) , T = T 

e e 


(5b) 

+ 

i, = 1 + e , 
1 m’ 

Pr + + 

b_ = 1 + “ — c , e = 

2 Pr^. m m 

e /» 
m 

(6) 


2.1 Turbulence Model 

The presence of e in b, and b„ and Pr. in b„ requires a turbulence model and 
in l c x c 

various expressions can be used for this purpose. The choice lies between an 
algebraic formulation in which detailed flow features, including transition 
from laminar and turbulent flow, are represented seml-empirically and models 
where the transport of turbulence properties is facilitated by differential 
equations. The models of the latter type, including the ubiquitous k-e 
model, offer the possibility of some degree of generality but, so far, have 
been unable to deal with low Reynolds number phenomena. Since the representa- 
tion of transition is essential to the present problem, we have preferred to 
use the algebraic eddy-viscosity approach of Cebecl and Smith, which is known 
to represent transition adequately. The formulation is given by a two-layer 
model defined by 
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(7a) 


{0.4y[l - exp(-y/A)]} 


3u . 

3y y tr - y - y c 


m 


0.0168 


J (u fi - u)dy 


Y tr y 


y c 5 y < 4 


(7b) 


where 


A = 26vu , 

T 


. w.1/2 
u = ( — ) 

T p 


Y * 


1 


1 + 5 . 5 ( y/6 ) 


(8) 


The condition used to define y Is the continuity of the eddy viscosity; from 

c 

the wall outward, Eq. (7a) Is applied until Its value Is equal to that given 
by Eq. (7b). 


In Eq. (7), y. Is an Intermlttency factor which accounts for the transl- 
t r 

tlonal region that exists between a laminar and turbulent flow. It Is given by 


Y tr = 1 - exp[-G( x - x tr ) J 4A ] 

x tr e 


(9) 


Here x^ r Is the location of the start of transition and the empirical factor 
G Is 


G = 


3 

J je p-1.34 

1200 2 x. 

v tr 


where the transition Reynolds number 



(u e x/v) 


tr’ 


( 10 ) 


Equation (9) has been determined In terms of flows where the freestream turbu- 
lence was characterized by low Intensities and small scale. Its extension to 
turbine blade passages which are subject to onset flows with turbulence Inten- 
sities In excess of 15% and scales larger than the width of the blade passage 
Imply that modifications may be required. The flow of Ref. [20] Involves a 
freestream Intensity In excess of 1% and Intensities of 10% In the near wakes 
of the circular cylinder which were used to simulate blade passing flows. As 
a consequence of the above, the Influence of Y^ r must be determined empir- 
ically and with consideration of these high turbulence flows. 
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2.2 Initial and Upstream Conditions 


The determination of the Initial conditions required for the above system Is 
Important and sometimes can be arbitrary but In that event, the values of 3u/3t 
at t = 0 are nonzero; this Implies an 1 nvl sc 1 d acceleration and, as a conse- 
quence, a slip velocity develops at the wall and Is smoothed by an Inner 

1/2 

boundary layer Initially of thickness (\>t) In which viscous forces are 
Important. Thus a double structure develops In the boundary layer and may be 
treated by the numerical method described In [28]. However, If Interest Is 
centered on the solution at large times, this feature may be reduced In impor- 
tance by requiring that the Initial velocity distribution satisfies the steady- 
state equation with the Instantaneous external velocity. In addition, it Is 
necessary to smooth the external velocity u g (x,t) so that 3u g /3t = 0 at t = 0 
and standard numerical methods may be used and remain stable. The use of a 
smoothing function Introduces some loss of accuracy at small values of t but 
the error soon decays to zero as the required value of u g Is approached. 

The calculation of upstream boundary conditions In the (t,y) plane at some 
x = x q when the conditions at a previous time line are known, can Introduce 
different problems. To Illustrate these difficulties for the case of a moving 
stagnation point, let us consider Eq. (1). Since u g * 0 at the stagnation 
point by definition. Its location, x $ , based on the external streamlines Is 
given by 


x $ = A(x $ ,t) (16) 

Figure 3 shows the variation of the stagnation point with time according to Eq. 

(1) with B(t) = 1 + c slnwt, c = 1 , u> = ir/4. We see that the stagnation point 

x Is at 2 when t = 2 and at 0 when t = 6, etc. If x were fixed, we could 
s 's’ 

assume that u = 0 at x $ = 1 for all time and for all n, but this Is not the 
case. It Is also possible to assume that the stagnation point Is coincident 
with zero u-veloclty for a prescribed time but we should note that the stagna- 
tion point given by Eq. (16) Is based on vanishing external velocity. For a 
time-dependent flow, as we shall discuss later, this does not necessarily 
Imply that the u-veloclty Is zero across the layer for a given E-location and 
specified time; flow reversals do occur due to the movement of the stagnation 


7 



-2 -1 Xs 0 1 

Figure 3. Variation of the stagnation point with time according to Eq. (1). 

point and cause the locus of zero u-veloclty to vary with x requiring the use 
of a special numerical method, as will be discussed In Section 3. 


2.3 Transformed Equations 

It Is more convenient and useful to express Eqs. (2)-(5) In a form more suit- 
able for computation. To achieve this, we define the dimensionless variables 
t, n, w, m, G together with a dimensionless stream function f. 


T 




u 1/2 

<*> > 


w = 


m 



d < T w - v 
y d 5 


Equations (2) to (5) become 



( 11 ) 

(u o vL) 1 / 2 f(^,n,T) 


where 


(b^f ") ' + f "8 + 

3w 3f 1 r , 
W 3E" 3 t + f 

3f 1 
35 

(12) 

(b 2 G') ' + G'0 + m( 1 

3G 

g) f 1 

36 

35 

(13) 

II 

O 

-h 

= f' = o. 

G = 0 

(14a) 

n . n e . f 

' = w, 

G = 1 

(14b) 

0 

3f 
= 3$ 


(15) 


As we shall discuss In Section 3, the upstream conditions can be obtained by 
solving the above system with a special numerical solution provided that the 
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Initial conditions are known. The latter can be obtained by making use of the 
steady-state conditions and expressing Eqs. (12) and (13) In the form 


' 8w df 1 

(b-|f ") * f-e + w-= f — 

( b 2 G' ) ' + G'6 + m(l - G)f ' = f ' || 

There Is no problem with the Initial conditions for these equations since 
calculations are first performed for laminar flows which admit similarity 
solutions at the stagnation point. The transitional and/or turbulent flow 
calculations can then be started at any location away from the stagnation 
point by specifying the location of the onset of transition. 


(16) 

(17) 
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3.0 NUMERICAL METHOD 


There are several numerical methods that can be used to solve the equations of 
the previous section, including the finite-difference methods due to Crank- 
Nicolson [29] and Keller [30]. Either scheme can be used satisfactorily for 
time-dependent boundary layers when there is no flow reversal across the layer 
although the Keller scheme has more desirable features than that of Crank- 
Nicolson for turbulent flows. When reversal occurs, and for problems associ- 
ated with upstream boundary conditions such as those considered here, the 
advantages of the Keller scheme become more pronounced and necessitate its use 
as we shall discuss here. Only the solution of the continuity and momentum 
equations are discussed. The energy equation is linear and is uncoupled from 
the momentum equation so that its solution is straightforward and Is described 
in detail in [27 ] . 

In this report the box scheme is used in two forms which depend upon the nature 
of the flow and, in particular on the presence or otherwise of reverse flow. 

The regular box scheme is used where the flow Is always in the streamwise 
direction and is described in the following subsection. The characteristic box 
scheme, developed by Cebeci and Stewartson, is used to compute flows with nega- 
tive streamwise velocities and also to generate the upstream boundary condi- 
tions discussed in subsection 2.2. 

3.1 Regular Box Scheme 

To solve the momentum equation for steady-state conditions subject to its 
boundary conditions by Keller's scheme known as the Box method, we let f ' = e 
and Introduce a new function g defined by 

e 1 = g (18a) 

and write Eqs. (15), (16) and (17) with b^ = b as 

0' = ff (18b) 

( b-j g ) 1 + ge + w = e (18c) 

n = 0, e = 0 = 0; n = n e ’ e = w (19) 


10 



To write the difference equations for the system given by Eqs. (18) and (19), 
we consider a net rectangle in which the net points are denoted by 


5 o - °- ■ E 1 -, * v 


"o - °- v Vi * v 


r = 1 , 2, . . . , I 


j - 1 , 2 , ...» J 


( 20 ) 


The finite-difference approximations to Eq. (18a) are obtained by averaging 
about the midpoint (^ , , 


.- 1 , 1 1 , 1 _ 1/0# 1 1 . 
h j (e j " Vl* " 9 j-l/2 " 1/2(9 j + ^-1* 


( 21 ) 


and those to Eqs. (18b, c) by centering all quantities except 0 at the center 
of the rectangle (l^ ^ ^ /2 ) and taking the values of each say e, 
at the corners of the box, that is, 

1-1 


1 - 1/2 1,1 1-1 x 1,1 1 1-1 1 - 1 , 
6 j — 1/2 * 2 e j-l/2 + C j-l/2 = 4 (e j + e j-l + e j + e j-l } 

The centering of 0 is achieved as 

1-1/2 1 1-1/2 J-1/2 = 1 fl , 

0 j-l/2 = 2 (0 j 6 j-l ) ~ 2 (0 j 0 j-l } 


(22a) 


(22b) 


In this notation, the difference approximations to Eqs. (18b, c) can be written 
in the form: 


h j ' 8 j-l ) = r 1 (e j-l/2 e j-l/2 ) (23a) 

h j ( Vj ‘ b J-l 9 j-l ) + 9 j — 1 /2 9 j — 1/2 + 9 j-l/2 6 j-l/2 ' r 1 (e ^ j-1/2 = R j-1 /2 


where 


(23b) 


R 


1-1 


-1, ,1-1 1-1 J-l 1-1, 


"I r . 2 , 1 


2 , 1-1 


2 , 1-1 


j-1/2 = " {h j (b j 9 j " b j-l 9 j-l ^ " r 1 [(W ^ j-1/2 ' (W ^ j-1/2 + (e ^-1/2^ 

(24) 

The boundary conditions, Eq. (19), become: 


e o = ®o - 


e J - W J 


(25) 


The system given by Eqs. (21), (23a, b) and (25) can be linearized by Newton's 
method and can be written as 




4e J - <e j-l - 2 (SS J * ’ <r l>j 


(26a) 
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(26b) 


40 J ' <9 J-1 - 2 <^>< Se J * ■ (r 3>M 

<S 1>J 49 J * (S 2>j Sg M * <s 3>j ir j " (S 4 ) j Sr 3-l * (s 5>J <e ] 

* (s 6>j‘ e M ' (l Vj <26c) 

6e 0 = 66 0 = 0, 6ej = 0 (27) 

where 


(r l } j = - (e j - Vl } + h j g J-l/2 (28a) 

(f 2 } j = R j -1/2 ' [h j (b j 9 j " b j-l g j-l ) + 2g j-l/2 9 j-l/2 " r 1 (e } j-l/2 ] (28b) 

(r Vj-l = " (6 j ' 6 j-7 ) * rj (e j-1/2 _ e j-1/2 ) {28c) 


(S 1 ) J = b j h j + 2 9 j -1/2 

(S 3>J 

1-1/2 
= 9 j -1/2 

■ 

(S 2 } j = ~ b j-l h j * 2 9 j — 1/2 

(S 4>3 

- {S 3>j 

<Vj • 


The resulting linear system, Eqs. (26) and (27), can be written In the form 



(30) 


and can be solved by the 3x3 algorithm based on the block elimination method 
described In [27]. 


3.2 Characteristic Box Scheme 

The procedure of the previous sub-section provides a solution of Eq. (16), 
starting at the stagnation point \ = l- , for a prescribed velocity distribution 
w(£) at t(0) = 0 and, as a result, the solutions on the next time "line" can be 
determined by an explicit method. If we wish to avoid stability problems, how- 
ever, an Implicit method Is required and the generation of upstream boundary 
conditions for the momentum equation (and the energy equation) becomes more 
difficult. An accurate and efficient method for this purpose Is the character- 
istic box scheme which Is able to deal with the moving stagnation point and the 
flow reversals which result. 
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To describe the characteristic box scheme, we note the definition of local 
streamlines and write 


di 


di 

e 


( 31 ) 


If we designate distance In this direction by s and the angle that It makes 
with the T-ax1s by P, (see Fig. 4) then Eq. (18c), which In terms of new 
variables e and g. Is 


(bg) ' 


. 3w dw 3e 3e 


(32) 


can be written as 

(bg ) 1 + 96 + 37 + w ^ \ (33) 

where 

\ = / l + e 2 (34a) 

p = tan’^e (34b) 



Figure 4. Notation and grid for the Characteristic Box scheme. For details, 
see [25]. 

The solution procedure for the system given by Eqs. (15a, b), (16) and (33) Is 
similar to that of the previous sub-section with small differences due to the 
manner In which the difference equations are written for Eq. (33). Also, since 
now we are considering a three-dimensional form of the equation, we also con- 
sider the net points 

t= 0, i = i , + k n = 1, 2 N (35) 

o n n -1 n 
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in addition to those given by Eq. (20). The difference approximations to Eqs. 
(18a, b) are Identical to those given by Eq. (21), (23a) except that the super- 
scripts 1 on each variable also Include n. To obtain the solutions on the 
first time line i , let us assume that the stagnation line based on the 
external streamline is given, the solution of point 1 Is known, and direct our 
attention to point 2. The finite-difference approximations to Eq. (33) are 
written along the streamline direction (see Fig. 4) at point P 


-1 


-1 


l 1 ,n1 ,rt\ .1 /t,m,n-l m,n-l 
(b j 9 j " b j-l g j-l J + o ( b i 9i 


2 '-j 

1 /J." 


h m » n - 1 m, n - 1 , 

Vl 9 M ’ 


dw 


dw 


+ 2 (9 j-l/2 + 9 j-l/2 )0 j-l/2 + 3 t (P) + W 3E (P) 


1/J.n ,m,n-l . "j-1/2 Hz 

2 j - 1 / 2 J-1/2 As j_l/ 2 


1,n m,n-l 
e - p 

^ i 1 /O C 1 1 


where 


(36) 


As 


j-1/2 ~ c °s3j_ 1/2 


(37) 


The relation between 0^ -j and those values of 0 centered at (1-1/2, n-1/2) 
and ( i-3/2, n-1/2) are J 


1-3/2 1-1/2 

0 P _ j-1/2 . ~ j-1/2 / rP c w gl-3/2 
J-1/2 ‘ 5, 3/2 - 5, 1/9 U " *1-3/2' j-1/2 


■ 1 - 1/2 


(38) 


m,n-l 


Rearranging Eq. (36) and denoting R^’.^ by 


r: 


m,n-l 


-l..m,n-l m,n-l ,m,n-l m,n-l 


3w 


3w v 


j-1/2 ■- [h J <b J 9j -bj_, gj _, ) * 2 aV <P) » 2(w ajXP) 

m,n-l 

,m,n-l e j-l/2 1 

Vl/2 A Sj _, /2 ] 

It can be shown that the linearized version of the resulting expression can be 
expressed In the same form as Eq. (26c), with 


(39) 


( S 1 ) j 


.-1 

h j C 7 


(S 2 ) j = _h j * C 7 


(40) 


/ v / r v 1 / i f n n) , n*-l i 

(S 3>j - (s 4 } j " ' 2 (9 j-l/2 9 j-l/2 ) 
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(s 5>j = - c 4 + 

J-l/2 


(s 6 } j = ~ c 5 + 1,n e j-l 

j — 1/2 


and 


. . n m,n-1 r , -1., 1,n 1,n ,1,n 1,n. . i,n m,n-l.-P 

r 2 j = R j -1/2 ‘ [h j (b j 9 j " VlVl* + (g j-l/2 + g j-l/2 )0 j-l/2 


1,n m,n-l 

,,1.n . e 1-)/2 .1,n e 1-l/2 ■ 

" < X J-V2 Vl/2> - x J-,/2 


where C] to c-j are given by 


E P - E 


c, = 


1-3/2 


1 ^ 1 - 3/2 " **1-1/2 


r J-3/2 
C 2 = r j-l/2 ’ 


C 3 ' As 


j - 1/2 


1 ,n 

c„ e 

C 4 = ~2 x 

j-l/2 


V,n e j — 1/2 j — 1 /2 j-l/2 J ’ 

- 1/2 

1 ,n 

r _ —2^ r — ini T,m»n-1 , 

c 5 ‘ 2 kl,n e j -1/2 j-l/2 j-l/2-' ’ 

j-l/2 


(41) 


(42) 


c = — e m,n_ ^ c - - fc (c - 0^ - ^) + c 1 
C 6 2 e j -1/2 ’ C 7 2 lC l 1C 2 j-l/2' C 2 J 


The resulting equation, which Is analogous to Eq. (23c) can then be solved 
together with Eqs. (23a, b) and (24) by using the same block-elimination method 
employed to solve the system given by Eqs. (23) and (24). 


To generate the upstream boundary conditions, let us Initially assume that the 
first profile on either side of the "edge" stagnation point Is known and that 
we are Interested In computing the flow say at 2 (see Fig. 4). If there Is now 
flow reversal at that point, then an extension of the procedure discussed In 
the previous sub-section can be used to obtain a solution there. In the pres- 
ence of flow reversal, we use the characteristic box scheme on an Iterative 
basis. To explain this further, let us assume that 6^ In Eq. (36) Is 
known, then the solution of the momentum equation at point 1 by the character- 
istic method can be achieved by solving Eqs. (23a, c) In the form 
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(43a) 


Se j ' Se j-1 " 2* (Sg 3 * <g 3-l ) " (r l*j 

(S 1>J 49 J * (s 2 > 3 i9 j-l f (s 5>J Se j * * (r 2> 3 (43b) 

subject to the boundary conditions 

6e 0 = 0, 6ej = 0 (44) 

Once a solution for the point 1 Is known, solutions to all the points to the 
right and to the left of point 1 can be obtained by the general characteristic 
scheme without further approximations. Thus the only approximation on a given 
time line results from the solutions generated at point 1. Our calculations 
have shown that the Influence of point 1 on its neighbors decreases with dis- 
tance. Thus the approximate solution at point 1 can be further Improved by 
repeating the calculations at point 1. 
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4.0 RESULTS AND DISCUSSION 


The numerical procedures of the previous section have been used to obtain solu- 
tions of the continuity, momentum and energy equation for three model problems. 
The results are presented and discussed In the following three subsections and 
correspond to uniform wall temperature (m = 0) and Prandtl number 0.72. The 
first two model problems Involve laminar flows and the third addresses laminar 
and turbulent flows. It should be emphasized that these model problems have 
been designed to test the numerical procedure which has been formulated In a 
general way so that Its application can readily be extended to the more prac- 
tical flows discussed In the Introduction. 

In regions where there was no flow reversal, the standard box scheme has been 
used. Where the calculations revealed flow reversal, the solution procedure 
subsequently made use of the characteristic box scheme. 

4.1 Oscillating Airfoil Model 

The first model problem corresponds to flow over an ellipse with a thickness 
ratio y(Hb/a) and with y « 1 at an angle of a. The surface of the body Is 
defined by 


x = -a cos<t>, y = ay s1n<t>, -ir < <J> < ir 


(45) 


With these definitions and to a first-order approximation, the external veloc- 
ity for the steady flow In the leading of a thin ellipse can be deduced from 
Invlscld flow theory to be 


u e (z) 


Z + Z 

0 


A77 


M6) 


Here u e (z) denotes a dimensionless velocity, the parameter z denotes a dimen- 
sionless distance related to the x- and y-coordlnates of the ellipse by x + a 
2 2 ” 2 

= (l/2)ay z , y = ay z measured from the nose, and z Q represents a reduced 

angle of attack. The parameter z Is also related to the dimensionless surface 

2 z 2 1/2 

distance \ with L = y a by E = J Q ( 1 + z ) dz. In the present study we extend 
Eq. (46) to unsteady flows by Introducing time dependency as 


u e (z * T) 


z + z q (1 + c slnwi) 


/7T7 


(47) 
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Figures 5 to 8 show the calculated results for two circular frequencies, u> = 
ir/30, */3 with z Q = 1 and c = -1/2. They demonstrate the effects of the 
frequency on the wall shear f"(0), wall heat flux G ' ( 0) and dimensionless 
displacement thickness, A*, which are defined by 




f"(0) = -5% /T , G ' ( 0 ) = - 

pu o k ^ T w-V 


4 * = uu s '^' n-S-ldy 

o 0 e 


(48a) 


(48b) 


The results of Figures 5a and 5b Indicate that the time dependent effects on 
the solutions Increase significantly as the circular frequency u> changes 
from u> =ir/30 to ir/3. While the wall shear and displacement thickness values 
computed at wt = ir and 2ir with w = ir/30 are nearly the same as those at t = 0, 
they differ considerably from each other when u = ir/3. The results In Figure 6 
show, however, that increasing the circular frequency by a factor of ten has 
practically no effect on the wall heat flux. 




Figure 5. The effect of frequency on the (a) wall shear parameter, f • ( 0) , (b) 
displacement thlcknesds. A*, of the first model problem. 
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Figure 6. Variation of the wall heat flux parameter, G'(0), with z for the 
first model problem with w = ir/30. The results with w = ir/3 
are virtually identical to those obtained with the smaller 
frequency. 



Figure 7. Variation of the velocity profiles for u> = ir/30 near the stagnation 
region of the first model problem for different values of wt, (a) 90, 
(b) 180, (c) 270, (d) 360. 
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Figure 8. Variation of the velocity profiles for u> = ir/3 near the stagna- 
tion region of the first model problem for different values of ojt, 
(a) 90, (b) 180, (c) 270, (d) 360. The dashed line Indicates the 
locus of zero u-veloclty. 


Figures 7 and 8 allow the examination of the effect the frequency has on the 
calculated velocity profiles In the vicinity of the stagnation point. Figure 
7 shows that the locus of the u-veloclty on time lines x = ir/2 and ir Is 
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essentially the same as In the steady case and as a result there are no flow 
reversals In the velocity profiles. However, as can be seen from Figure 8, 
Increasing the frequency to ir/3, flow reversals begin to occur around the 
stagnation point and become rather prolonged as time Increases to wt = 3ir/2. 
At wt = 2ir, the region of flow reversal Is reduced but Is not zero, as It was 
at wt =0. 


4.2 Gust Response Model 

The second model problem allows us to study the transient motion between the 
two steady stagnation flows similar to those an airfoil would encounter during 
a gust or change In Incidence as an Initially steady-flow accelerates from 
u /u = % to u /u = 2(5 - 1). For this purpose we use Eq. (1) with w = ir/8 

0 oo e a> 

and 


A(t) = 2 - cos^wt, B(t) = sln^wt 0 < t < 8 (49) 

This flow has been studied by Cebecl, Stewartson and Williams [28] with 
slightly different forms of A(t) and B(t) 

A(t) = 2 - e-\ B(t) = 1 - e~ T (50) 

than those used here. Our modified expressions were required to ensure that 
(du g /3t) Is equal to zero at t = 0, as discussed previously. 

Figure 9 shows a comparison between the two external velocities considered In 
[28] and In the present study, and Figure 10 shows the results of the present 



Figure 9. Variation of external velocity with time for the second model 
problem. 
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Figure 10. Variation with time of (a) wall shear parameter, f " ( 0 ) , (b) wall 
heat flux parameter, G'(0), and (c) displacement thickness, A* 
for the second model problem. 

calculations together with those of [28]. Our calculated results correspond 

to nondlmenslonal forms of the wall shear, f"(0), wall heat flux, G'(0), and 

★ 

displacement thickness, A*(5 /R^ 6 q /L) where 

1 °° 

K = n * ( u p " u)dy (51) 

0 u o 0 e 

The results of [28] were obtained for momentum based properties which are 
linked to the above definitions by 

f " ( 0 ) = (l - B)t + g"(0) (52) 

A* = (1 - B)6* - g a (53) 

★ 

with t, g"(0), 6 and g denoting here the parameters computed in [28]. 
o ® 

As expected, the results In Figures 10a and 10c show that the behavior of f"(0) 
and A* is different during the transient motion due to the use of different 
functions for A(t) and B(t). In each case, however, the solutions approach 
the asymptotic values relatively fast. The results for G'(0), shown In Fig. 
10b, were obtained with the present method since the study of [28] did not 
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consider the energy equation. Again, as expected, the solutions of the energy 
equation approach the asymptotic value of 0.7095 [27]. 

4.3 Rotor Wake Model 

The third model problem Involves flow on a circular cylinder of diameter D 
which experiences the periodic passing of wakes. It corresponds to the experi- 
mental arrangement of [20], see Fig. 2, and permits comparison between calcu- 
lations performed for turbulent flows with measurements. The model problem Is 
first addressed In terms of laminar flow to test the numerical procedure. 

Figure 11 shows the flow configuration and the notation for this model problem. 

For a time period t , the cylinder is subjected to a freestream velocity V and 

9 00 

for t It Is immersed in a superimposed rotating wake which has a rotational 

component wr. The cycle repeats Itself with a blade passing frequency F(El/(t^ 

+ t )] and Is related to the Strouhal number St by FD/V . To derive dlmenslon- 
w 00 

less expressions for the parameters A(£,t) and B(i) In Eq. (1) needed to define 

the external velocity distribution near the stagnation region, we adopt the 

procedure discussed In the Appendix. 



Figure 11. Notation and flow configuration for the rotor-wake model. 

At first the calculations were performed for a laminar flow In order to test 

the numerical procedure. With the choice of E„ = /TO, E, = 100 and V = 1/3, 

2 1 m 

the parameters A, B and f, Eqs. (A-15) to (A-17) become 

A(t) = [f 2 + 10(St) 2 (l - f) 2 ] 1/2 (54) 
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(55) 


B( t ) = tan 1 C v50 St 

f = | + 1 cos [ir( 1 - 100 Stx) ] 

The computed values of wall heat flux, G ' ( 0) for two values of Strouhal number 
show that they are not Influenced by the changes in the freestream velocity and 
are virtually constant for the range of £ and t values considered with G'(0) 

~ 0.50 for St = 0.1 and G'(0) ~ 0.51 for St = 0.2. On the other hand, as 
shown in Fig. 13, the computed values of wall shear, f"(0), for St = 0.1 are 
significantly Influenced by the changes in the freestream velocity which causes 
flow reversals in the velocity profiles around the stagnation point based on 
the vanishing of the external velocity. The movement of the stagnation point 
and the resulting flow reversals Increase with time and with space. For 
example, the calculations for steady state have the stagnation point at £ = 0, 
and, as expected, there is no flow reversal on either side of the stagnation 
point. At t = 0.05, the stagnation point moves to l = 0.15 but the flow rever- 
sals in the velocity profiles continue up to and including \ = 0.85 as can be 
seen from the results shown in Fig. 14. At t = 0.10, Fig. 13, the stagnation 
point has moved to \ = 0.55 but the flow reversals persist for a longer distance 
and continue until £ = 1.20. As can be seen from the velocity profiles in Fig. 



Figure 13. Variation of the wall shear parameter, f"(0), with \ for the 
third' model problem. St = 0.10. 
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Figure 14. Variation of the velocity profiles near the stagnation region of 
the third model problem for two values of time, (a) t = 0.05, 
(b) t = 0.10. 


14, the region of flow reversal across the layer has now Increased and Is sub- 
stantially more pronounced than that at t = 0.05. 


The wall shear and displacement thickness results for St = 0.2 show similar 
trends. As expected, the flow reversals In the velocity profiles, for example 
at t = 0.10, are bigger than those at St = 0.1 but cover the same range In £. 

The above laminar flow results show that the numerical procedure Is able to 
obtain solutions for a range of blade-passing frequencies of practical rele- 
vance. The movement of the stagnation point with space and time and the 

resulting flow reversals around the stagnation point cause no computational 

l 

difficulties and the numerical tests show that the accuracy Is better than 
required for practical problems. 

The present method Is applicable to laminar, transitional and turbulent flows. 
It requires that the calculations start as laminar from the stagnation point. 
The onset of transition can be assigned to occur In accordance with experiments 
and the subsequent transitional and turbulent flows can be represented by the 
algebraic eddy-viscosity model described In Section 2.2. As a consequence of 
this model, the length of the transitional region Is depended upon the free- 
! stream conditions. 

i In the case of turbulent flow, estimates of E., , E„ and V can be obtained from 

12 m 

the experiments of Ref. 20. We note that D represents the diameter of the 
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heated stator, d the diameter of the rotating pin, S the separation distance 

between the two pins of diameter d and D (see Fig. 2), r the distance from the 

axis of rotation, C Q the drag coefficient of the n number of rotating pins, 

and V the wake centerline velocity. The experimental values are D = 12.7mm, 
m 

d = 3.18mm, S = 25.4mm, r = 16.925cm, C = 1.2, n = 24, V = 0.62, so that 

D m 

E ] = 79, E ? = 3.49 and 


2 2 2 2 1/2 

A(T) = [r ♦ (3.49) (St) (1 - r)Y ,d (57) 

B(t) = tan" 1 [3.49 (^-j-^)] (58) 


f = 0.81 + 0.19 cos[ir(l - 79 Stx)] 


(59) 


Calculations using Eqs. (57) to (59) In Eq. (1) were performed In the above 
manner to simulate the near stagnation region of the flow of [20], The onset 
of transition was set very close to the stagnation point and transitional and 

turbulent flow calculations were performed for 0 < t < t + tg, for one cycle 

w 

for the experimental Reynolds number of R. = 76,000. Since the experimental 

d 

data suggests that the average Nusselt number Is relatively constant and Is 
ten-percent higher than that of laminar flow, It was assumed that the transi- 
tional region was negligibly short. Furthermore, the present eddy-viscosity 
model Is a function of the velocity field and the heat transfer parameter Is 
relatively constant In the stagnation region 0 < 9 < 40°, so that It was fur- 
ther assumed that the eddy viscosity In this region be constant. This implies 
that the distribution of eddy viscosity determined at the end of the transi- 
tional region [y^ r = 1] retains the same numerical value throughout the turbu- 
lent flow calculations. As a consequence of the above, the ratios of the 
aver- age Nusselt number, Nu, on the blade passing a wake to that of a blade 
In a freestream during one cycle, Nu $ , used In the presentation of the 
experl- mental results were calculated from 



Nu k + Nu„ 

b a _ 

Nu. 


t 

J w Nudt + 


t +t 


w 


Nu^dt 


w 


t +t 

J ® w Nu„dt 
0 1 


(60) 


Since the cycle time Is t^ + t , the numerator of Eq. (60) contains two parts, 

one for the blade being submerged In the wake, Nu b , and the other for the base 

In a freestream, Nu , during time t where the flow Is laminar and admits 

9 9 
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similarity. In the transformed coordinates of the calculations, the ratio of the 
averaged Nusselt number may be written as 


Nu 

Nu^ 


/V * <W 9 

<W‘g * V 


(61) 


I 

where (G ) 0 = 1.0034 for laminar flows. The measured and calculated values of 
average Nusselt number ratios were found to be 1.1 and 1.09, respectively, and 
constant throughout the turbulent flow region. The closeness of the two 
results Is gratifying but should be viewed with caution bearing In mind the 


assumptions which have been made In relation to the eddy-viscosity model and 
the assumed freestream velocity distribution. Alternative approaches to the 
representation of the turbulence characteristics of the stagnation region have 
been examined, for example by Taulbee and Tran [31], and deserve future con- 

I 

slderatlon. In addition, the integration of the heat transfer parameter G w 
Is subject to some uncertainty In the blending region between laminar and 

turbulent flows. 
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5.0 CONCLUDING REMARKS 


The work of the preceding sections represents essential steps In the develop- 
ment of a general method for the calculation of unsteady heat transfer on 
turbine blades. It has emphasized the stagnation region since this must be 
correctly represented In order for the method to calculate the flow and heat- 
transfer characteristics of a blade passage. This region Involves a moving 
stagnation point with consequent reverse flows and has required the development 
and use of a novel numerical procedure to solve the equations for conservation 
of mass, momentum and energy. The characteristic box scheme, with Its stabil- 
ity requirements, has been used In regions of reverse flow and the standard box 
scheme elsewhere. Transitional and turbulent flow have been represented by an 
eddy-viscosity formula. 

The method has been applied to three model problems devised to allow Its quan- 
titative evaluation. They correspond to an ellipse, with leading-edge geometry 
similar to that of a thin airfoil, oscillating with a uniform-velocity onset 
flow; a stationary airfoil subjected to an onset velocity which changes sud- 
denly; and to a cylinder which experiences an onset flow with periodic wakes. 
Laminar-flow calculations were performed for these three model flows and 
results presented to demonstrate that the method is able to calculate their 
essential features In a convenient and numerically accurate manner. The third 
problem simulates an experimental arrangement for which heat-transfer measure- 
ments have been reported and calculations were performed to Include the transi- 
tional- and turbulent-flow regions; again the method predicted results which 
represented the essential features of the flow. Including an averaged Nusselt 
number in close agreement with the measured value. 

The results confirm that the calculation method correctly predicts the flow 
characteristics of the stagnation region of unsteady flows of relevance to 
blades subjected to onset velocities which vary In space and In time. The 
method can now be extended to Include procedures for the solution of the 
unsteady, Invl scld-f low equations and to ensure Interaction between the Invls- 
c 1 d and boundary-layer flows. Both procedures have already been developed so 
that the general method can be assembled and applied to blade-passage flows. 
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Appendix 

EXTERNAL VELOCITY DISTRIBUTION IN THE STAGNATION REGION 
WITH NONUNIFORM ONSET FLOW 


Consider a fixed circular cylinder of diameter D In a uniform flowfleld V 

oo 

which Is disturbed by passage of wakes shed from rotating blades with a rota- 
tional velocity wr. Assume the wakes to be nonoverlapping so that a single 
wake (for simplicity taken as symmetric) relative to the blade Is given as 
shown In Fig. 11. To determine the onset flow V and the angle of attack a 
with respect of the stationary cylinder as a function of time, we observe that 
the onset flow velocity relative to the moving cylinder Is and that the 
cylinder sheds a wake whose velocity distribution parallel to the V -dlrec- 

K 

tlon Is given by 

v i ■ V<«> 


Here y denotes the distance normal to V and 6 the wake half-width. The 

K 

velocity In the wake with respect to the stationary cylinder, V, can then be 
obtained by subtracting the rotational speed wr from V 1 as shown In Fig. A1 . 
In terms of known quantities, It follows from Fig. A1 that 



2 2 
ii> r 


(1 


f)2 ] l / 2 


( A-l ) 


4- t * >r 

tana = — 


00 



( A- 2 ) 



Figure A1 . Sketch of velocity vectors behind a moving blade row. 

To apply Eqs. (A-l) and (A-2) to the stagnation point solution on the circular 
cylinder, we make further simplifying assumptions as follows: although Eqs. 
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(A-l) and (A-2) represent a variable onset flow for the whole cylinder, we take 
the stagnation point location to be determined by the angle a, or that 0 = a 
on the cylinder periphery, and let the Instantaneous onset flow be equal to V. 
This Is clearly a quasi-steady approach which becomes more accurate as the 
ratio of the cylinder diameter with respect to the wake width decreases. By 
analogy to the steady-state flow around the stagnation point, we write for the 
velocity at each Instant 

“e ■ D* (X - V (A - 3 > 

where X Is the stagnation point location with respect to a = 0 = 0. 
o 

Inserting Eqs. (A-l) and (A-2) Into Eq. (A-3), we obtain: 


u .22 

f = 2[f 2 ♦ (1 

V CO V 


f> 2 ] 1/2 ( *§ 


tan' 


r— (1 
L V 1 


-)]} 


( A-4) 


For a symmetrical wake, the nondlmenslonal velocity distribution f may be 
represented by: 

C d 1/2 .3/2 2 

f({) = 1 - 0.976 (-jj-) [1 - ||| ] ( A- 5 ) 

Here C^d Is the drag area of the blade and S Is the distance from the blade 
to the cylinder. Equation (A-5) can also be approximated by 


1 + V 


f (f) - 


m 


1 - V 

- - - cosir (*) 


with V m denoting the nondlmenslonal centerline velocity given by 


V 1/2 

v m ■ 1 - °- 976 <-T> 


To Introduce time dependence, we observe that the wake width W In the plane of 
rotation of the cylinder Is: 


. 26 _ 

slny 



(A-6) 


Since one wake covers 6 V^/irr of the periphery, the period t^ for 
one wake to pass Is: 

2 2( [1 * Cr/VJ 2 ] 172 
t w ~ rwV = V (wr/V ) 

CO 00 00 ' 


(A-7) 
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If there are n blades, the period tg for one gap Is 

t .fc.t 

g nu w 

or In terms of the blade passing frequency F 

nw 


F = 


2ir 


The full width of the wake, 26, Is given by 


26 = 1 . 1 32( CpdS) 


1/2 


It takes time t to travel across the wake of distance 26, 

w 

t, the normal distance traveled Is: 


- 26 1 


y = 


w 


Here y Is measured from the edge of the wake 


y * 5 - y 


so that: 


Jt _ 1 _ Jt _ 1 _ 2£ 
4 4 *w 


0 £ t < t w 


In terms of the relations defined by Eqs. (A-7), (A-9), (/ 
with the definition of Strouhal number St, Eq. (A-12) can 


* . t . E 
6 ' t l 


where t Is a dimensionless time 


T St 


[1 + E 2 2 St 2 ] 1/2 


T = 


V t 

oo 

IT 


and the parameters E-| and E 2 are defined by 

2ir 


E 1 - 


. 1 32n( C D d/0 S/0) 


1/2 


2irr 
'2 = nD 


( A-8a) 
( A- 8b ) 
(A-9) 

so that In time 

(A-10) 

(A-ll) 

(A-12) 

-10) and (A-ll), and 
also be written as 

( A— 1 3 ) 

( A-l 4a ) 
(A-14b) 
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Noting that f Is a function of i and that Eq. (A-4) is of the form given by 
Eq. (1), we can represent the external velocity distribution near the stagna- 
tion region of the circular cylinder with u = 2V , \ = 2x/D = x/r, St = (n/ir)(u>r/ 

00 00 

V ) by defining A and B as 
00 

2 2 2 21/2 

A(5.t) - [f + E ^t) 0 - f) ] ( A— 1 5 ) 

B(i) = tan -1 [E 2 St ( I ~)] ( A- 1 6 ) 

1 + V 1 - V 

f(T) = 2 ~ + — 2 “^ cos [ir (y/4) ] (A-17) 

with y/6 In Eq. (A-17) given by Eq. (A-13). 

As we can see from Eqs. (A-15), (A-16) and (A-17), an expression of the form 
given by Eq. (1) can be obtained for the stagnation region of the cylinder as 
a function of St by assigning values to E 1 , E,, and V . This can be done 
for both laminar and turbulent flows as discussed In Section 4.3. 
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